Quartically regularized polynomials
Adapted from: Section 3.3.1 of (Cartis and Zhu, 2026)
Higher-order optimization methods minimize, at each iteration, a cubic Taylor model regularized by a quartic term:
\[m(s) = f_0 + g^\top s + \frac{1}{2} H[s]^2 + \frac{1}{6} T[s]^3 + \frac{\sigma}{4} \|s\|_2^4.\]
For any $\sigma > 0$, this model is bounded below so $m(s) - m^*$ is nonnegative where $m^*$ is its global minimum. If $m(s) - m^*$ is moreover a sum of squares then $m^*$ can be computed exactly with a single Sum-of-Squares program, i.e., the zeroth level of the Lasserre hierarchy. (Cartis and Zhu, 2026) shows that this is the case for $\sigma$ large enough and that, when $T = 0$, it is the case for all $\sigma > 0$. They also show that the geometry of the regularization matters: with the separable quartic regularization $\|s\|_4^4 = s_1^4 + \cdots + s_n^4$ instead of the Euclidean one $\|s\|_2^4 = (s_1^2 + \cdots + s_n^2)^2$, $m(s) - m^*$ may fail to be a sum of squares for every $\sigma > 0$. In this tutorial, we verify this numerically on their example.
using LinearAlgebra
using DynamicPolynomials
@polyvar s[1:3](DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}[s₁, s₂, s₃],)We consider the following nonconvex quadratic part whose Hessian has eigenvalues $-4$, $-4$ and $20$:
quad = 2 * sum(s .^ 2) + 8 * (s[1] * s[2] + s[1] * s[3] + s[2] * s[3])\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} \]
We need to pick an SDP solver, see here for a list of the available choices.
using SumOfSquares
import Hypatia
import Dualization
solver = Dualization.dual_optimizer(
optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true),
)#7 (generic function with 1 method)The following function computes the largest $\gamma$ such that multiplier * (p - γ) is a sum of squares. With the default multiplier, this is the zeroth level of the Lasserre hierarchy. It returns the lower bound $\gamma$ and the constraint reference, which we will use to extract minimizers from the moment matrix.
function sos_lower_bound(p; multiplier = 1)
model = SOSModel(solver)
@variable(model, γ)
@objective(model, Max, γ)
con_ref = @constraint(model, multiplier * (p - γ) >= 0)
optimize!(model)
@assert termination_status(model) == MOI.OPTIMAL
return value(γ), con_ref
endsos_lower_bound (generic function with 1 method)Separable quartic regularization
We first consider the separable regularization with $\sigma = 4$. The resulting polynomial was shown in (Ahmadi et al., 2023; Theorem 3.3) to be such that $m_A(s) - m_A^*$ is not a sum of squares.
m_A = quad + sum(s .^ 4)\[ 2s_{3}^{2} + 8s_{2}s_{3} + 2s_{2}^{2} + 8s_{1}s_{3} + 8s_{1}s_{2} + 2s_{1}^{2} + s_{3}^{4} + s_{2}^{4} + s_{1}^{4} \]
The zeroth level of the Lasserre hierarchy gives the following lower bound:
γ_A, _ = sos_lower_bound(m_A)
γ_A-2.3859528350590975Is this lower bound tight ? Recall that a polynomial $p$ is nonnegative if $(s_1^2 + s_2^2 + s_3^2) p(s)$ is a sum of squares, see the Motzkin tutorial. Using this multiplier, we get a larger lower bound:
γ_A_mult, con_ref = sos_lower_bound(m_A, multiplier = sum(s .^ 2))
γ_A_mult-2.144278447395472We can extract the minimizers from the moment matrix:
atoms = atomic_measure(moment_matrix(con_ref), 1e-4)Atomic measure on the variables s[1], s[2], s[3] with 7 atoms:
at [-0.6674086334298758, 1.1497482276306912, -0.6674025653369381] with weight 0.07532816260714947
at [0.6674086667533291, 0.6674025714715247, -1.1497482877696992] with weight 0.0753038701744452
at [-1.1497467185458887, 0.6674041807780484, 0.6674042581873538] with weight 0.07532761207921562
at [1.1497466371206675, -0.6674041672390737, -0.6674041679530347] with weight 0.07532763395498196
at [0.6674070898031693, -1.1497476675493068, 0.6674030764542367] with weight 0.0753283128180053
at [-0.6674070417013989, -0.6674031450918838, 1.1497476864180787] with weight 0.07530403309330821
at [0.0, 0.0, 0.0] with weight 0.0We recover the six global minimizers $\pm(-1.1498, 0.6674, 0.6674)$ and their permutations. The moment matrix also contains an atom at the origin with zero weight. It comes from the multiplier which vanishes at the origin, so we ignore it.
minimizers = [atom.center for atom in atoms.atoms if atom.weight > 1e-4]6-element Vector{Vector{Float64}}:
[-0.6674086334298758, 1.1497482276306912, -0.6674025653369381]
[0.6674086667533291, 0.6674025714715247, -1.1497482877696992]
[-1.1497467185458887, 0.6674041807780484, 0.6674042581873538]
[1.1497466371206675, -0.6674041672390737, -0.6674041679530347]
[0.6674070898031693, -1.1497476675493068, 0.6674030764542367]
[-0.6674070417013989, -0.6674031450918838, 1.1497476864180787]Evaluating $m_A$ at the minimizers confirms that this is the global minimum and that the zeroth level of the hierarchy has a gap:
m_A_min = minimum(x -> m_A(s => x), minimizers)
m_A_min-2.1442784145215135As shown in (Cartis and Zhu, 2026; Lemma 3.3), the gap does not close as $\sigma$ increases. Indeed, with $\tilde{s} = \sqrt{\sigma}s/2$, we have $\mathrm{quad}(s) + \sigma \|s\|_4^4/4 = 4 m_A(\tilde{s})/\sigma$ so the ratio between the lower bound and the minimum is independent of $\sigma$:
for σ in [1, 10, 100]
m = quad + σ / 4 * sum(s .^ 4)
γ, _ = sos_lower_bound(m)
println("σ = $σ: lower bound = $γ, minimum = $(4m_A_min / σ)")
endσ = 1: lower bound = -9.54381088156368, minimum = -8.577113658086054
σ = 10: lower bound = -0.9543811375079798, minimum = -0.8577113658086054
σ = 100: lower bound = -0.09543811800874058, minimum = -0.08577113658086054Euclidean quartic regularization
Let us now replace the separable regularization with the Euclidean one:
σ = 4
m_E = quad + σ / 4 * sum(s .^ 2)^2\[ 2.0s_{3}^{2} + 8.0s_{2}s_{3} + 2.0s_{2}^{2} + 8.0s_{1}s_{3} + 8.0s_{1}s_{2} + 2.0s_{1}^{2} + s_{3}^{4} + 2.0s_{2}^{2}s_{3}^{2} + s_{2}^{4} + 2.0s_{1}^{2}s_{3}^{2} + 2.0s_{1}^{2}s_{2}^{2} + s_{1}^{4} \]
This is a quadratic model with Euclidean quartic regularization ($T = 0$) so $m_E(s) - m_E^*$ is a sum of squares for every $\sigma > 0$. Indeed, it can be checked that
\[m_E(s) + \frac{4}{\sigma} = \frac{\sigma}{4}\left(\|s\|_2^2 - \frac{4}{\sigma}\right)^2 + 4(s_1 + s_2 + s_3)^2.\]
The zeroth level of the hierarchy finds this lower bound $-4/\sigma = -1$:
γ_E, con_ref = sos_lower_bound(m_E)
γ_E-0.9999999899545772The global minimizers form the circle $s_1 + s_2 + s_3 = 0$, $\|s\|_2 = 2/\sqrt{\sigma}$. We use an approximation of rank 4 to extract a few points of this circle. This is a heuristic for finding minimizers; the resulting atomic measure need not reproduce all the moments:
atoms = atomic_measure(moment_matrix(con_ref), FixedRank(4))Atomic measure on the variables s[1], s[2], s[3] with 4 atoms:
at [-0.027343615746903256, -0.6930383487994694, 0.720381962467486] with weight 0.3003636768217744
at [0.02734361574690349, 0.6930383487994702, -0.720381962467487] with weight 0.30036367682177356
at [-0.7304699581580496, 0.6811625482830396, 0.04930740169524811] with weight 0.1939083213868171
at [0.7304699581580503, -0.6811625482830401, -0.0493074016952482] with weight 0.1939083213868167We can verify that they belong to the circle:
for atom in atoms.atoms
x = atom.center
println("sum = $(sum(x)), norm = $(norm(x))")
endsum = -2.0788866184062726e-9, norm = 0.9999999990386578
sum = 2.078886729428575e-9, norm = 0.9999999990386591
sum = -8.179761923354789e-9, norm = 0.9999999984084017
sum = 8.179761944171471e-9, norm = 0.9999999984084026Adding a cubic term
We now add a small cubic term $0.1(s_1^3 + s_2^3 + s_3^3)$ to both models. With the separable regularization, there is still a gap between the lower bound of the zeroth level and the global minimum:
cubic = 0.1 * sum(s .^ 3)
m_A2 = m_A + cubic
γ_A2, _ = sos_lower_bound(m_A2)
γ_A2-2.3989906372281125The multiplier closes the gap again and gives the global minimum:
γ_A2_mult, con_ref = sos_lower_bound(m_A2, multiplier = sum(s .^ 2))
γ_A2_mult-2.240948583236178With the Euclidean regularization, the model is now a cubic model so (Cartis and Zhu, 2026; Theorem 2.3) only guarantees that $m(s) - m^*$ is a sum of squares for $\sigma$ large enough. We compare the lower bound found by the zeroth level of the hierarchy with the value of $m$ at the minimizers extracted from the moment matrix. If they match, then the lower bound is the global minimum.
for σ in [0.1, 1, 4, 10, 40]
m = quad + cubic + σ / 4 * sum(s .^ 2)^2
γ, c = sos_lower_bound(m)
μ = atomic_measure(moment_matrix(c), 1e-4)
m_min = minimum(atom -> m(s => atom.center), μ.atoms)
println("σ = $σ: lower bound = $γ, value at extracted minimizers = $m_min")
endσ = 0.1: lower bound = -53.43764871528847, value at extracted minimizers = -53.437654123518456
σ = 1.0: lower bound = -4.353277454347697, value at extracted minimizers = -4.353277566345602
σ = 4.0: lower bound = -1.0424382420832403, value at extracted minimizers = -1.0424382409154571
σ = 10.0: lower bound = -0.410583054678908, value at extracted minimizers = -0.41058305172572457
σ = 40.0: lower bound = -0.10130678112733997, value at extracted minimizers = -0.10130677722939606For this example, the zeroth level is exact even for small values of $\sigma$. The theoretical threshold of (Cartis and Zhu, 2026; Theorem 2.3) is sufficient but not necessary.
This page was generated using Literate.jl.